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ABSTRACT 

Improving upon our purely dynamical work, we present three-dimensional simulations of the 
atmospheric circulation on Earth-like (exo)planets and hot Jupiters using the GFDL-Princeton 
Flexible Modeling System (FMS). As the first steps away from the dynamical benchmarks of 
Heng, Menou & Phillipps (2011), we add dual-band radiative transfer and dry convective 
adjustment schemes to our computational setup. Our treatment of radiative transfer assumes 
stellar irradiation to peak at a wavelength shorter than and distinct from that at which the 
exoplanet re-emits radiation ("shortwave" versus "longwave"), and also uses a two-stream 
approximation. Convection is mimicked by adjusting unstable lapse rates to the dry adiabat. 
The bottom of the atmosphere is bounded by a uniform slab with a finite thermal inertia. For 
our models of hot Jupiters, we include an analytical formalism for calculating temperature- 
pressure profiles, in radiative equilibrium, which accounts for the effect of collision-induced 
absorption via a single parameter We discuss our results within the context of: the predicted 
temperature-pressure profiles and the absence/presence of a temperature inversion; the pos- 
sible maintenance, via atmospheric circulation, of the putative high-altitude, shortwave ab- 
sorber expected to produce these inversions; the angular/temporal offset of the hot spot from 
the substellar point, its robustness to our ignorance of hyperviscosity and hence its utility in 
distinguishing between different hot Jovian atmospheres; and various zonal-mean flow quan- 
tities. Our work bridges the gap between three-dimensional simulations which are purely dy- 
namical and those which incorporate multi-band radiative transfer, thus contributing to the 
construction of a required hierarchy of three-dimensional theoretical models. 
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1 INTRODUCTION 

Much of the prev ious research on thre e-dimension al simulations of exoplanetary atmospheric circulation ha s focused on the atmo- 



spheric dynamics (e.g..lCooDer & Showmanlb oOS. 200^ 



Koskinen et alj|2007 



Showman et al 
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Dobbs-Dixon. Gumming & Li rjl20inl: 
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Selsis. Wordsworth & Forget.2011i) have developed sophisticated schemes which combine dynamics with multi-band radiative transfer. De- 



spite making progress, our understanding of atmospheric circulation on extrasolar planets (or "exoplanets") remains incomplete. To move 
forward, we need to construct a hierarchy of three-dimensional simulations o f varying sophistication, so as to e l ucidate the isolated e ffects 
of various pieces of physics as well as the complex interplay between them i Peixoto & Gort 1984 : Held 2005 : Pierrehumberll 2010[) . The 
purpose of the present study is to bridge the gap between these two bodies of work by elucidating the details of a computational setup of 
intermediate soph istication, so as to serve as a founda t ion for future w ork. Our work ma y be partially regarded as the three-dimensional 
generalizations of iHubenv, Burrows & SudarskvlfcoOSb . lHanserj J2008h and lGuilloj j2O10h . 
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As natural extensions to purely dynamical simulation work, the effects of radiative transfer and convection need to be considered. As a 
first step in adding radiative transfer, we simplify the treatment by assuming that most of the stellar emission peaks at a wavelength w hich 
is shorter than and distinct from that at which the exoplanet re-radiates the absorbed heat (Hansen 2008; Langton & Laughlin 2008; Guilloll 



20 1 d; IPierrehumbertlbo 1 ol) . For example, the Sun-like host star of an exoplanet has an emission spectrum which peaks in the optical, since 



A,_^0.5Mm (1) 

using Wien's law, where T* is the effective temperature of the star. For M stars (red dwarfs), the stellar spectrum instead peaks at A*, max > 
0.8 /^m (3700 K/T*). The equilibrium or effective temperature of the exoplanet i^ 

T.^ 1300Kr^U|lW^)-^^ (2) 



VeOQQKy' \Rq) V0.05AU/ 
where is the stellar radius and a denotes the spatial separation between the star and the exoplanet. As the stellar photons travel downward 
into the atmosphere, they are eventually absorbed and re-emitted. The emission spectrum of an exoplanet is expected to peak in the infrared, 
since 

Thus, to a first approximation one can regard the star as emitting "shortwave" radiation onto the exoplanet, which then re-emits in the 
infrared ("longwave"). The conversion of a shortwave photon into multiple longwave photons serves to increase entropy. Our treatment of 
the radiative transfer then requires the specification of the shortwave (rs) and longwave (tl) optical depths, as well as the flux associated 
with stellar irradiation ( ). Also as a fir st step, we mimic convection by adjusting unstable atmospheric lapse rates to the stable dry adiabat 
i Manabe. Smagorinskv & Strickleiill965l) . The key advantages of our simpler setup, compared to simulations with multi-band radiative 
transfer, are that it allows for a more efficient sampling of parameter space and clean comparisons to analytical models. 

In [J2I we describe our computational setup, including our schemes for stellar irradiation, shortwave/longwave absorption, radiative 
transfer and convection. Prior to simulating the atmospheric circulation on hot Jupiters, we test our computational setup using Earth-like 
models (©. In g] we describe our method for simulating hot Jovian atmospheres — including an analytical formalism for calculating 
temperature-pressure profiles which includes the effect of collision-induced absorption — and examine results from several models based on 
the values of parameters appropriate to HD 209458b. The implications of our results, as well as a concise summary of the present study, are 
discussed in !j5] Table [T] summarizes the list of adopted parameter values as well as the commonly-used symbols. Appendix lAl contains the 
technical details of the atmospheric boundary layer scheme used only for our Earth-like models. 



2 THE GFDL-PRINCETON FLEXIBLE MODELING SYSTEM 

In this section, we describe our computational setu p, which is based on the Lima release of the FMS, developed by the Geophysical Fluid 
Dynamics Laboratory (GFDL) at Princeton University JOordon & Stemiri982 ). As the FMS is set up to work in MKS (metres, kilogrammes 
and seconds) units, most of the discussion will follow suit. Unless specified otherwise, the term "day" refers to an Earth day (exactly 86400 
s). 



2.1 Dynamics 

We implement the spectral dynamical core, which solves the three-dimensional primitive equations of meteorology — these equations 
assume hydrostatic equilibrium and other consistent approximations such as a small aspect ratio for the atmosphere (e.g., see Heng. Menou & Phillippd 



20 111 and references therein). Held & Suarez (1994.) originally proposed the comparison of purely dynamical simul ations with simplified ther- 



mal fo rcing and performed via different methods of solutions, known as the "Held-Suarez benchmark" for Earth. iHeng. Menou & PhillippsI 



1 201 Ih extended the Held-Suarez benchmark to tidally-locked exoplanets by examining three additional cases: a hypothetical exo-Earth, a 



shallow model for hot Jupiters and a deep model of HD 2O9458b0 The spectral and finite difference (B-grid) dynamical cores of the FMS 
were both subjected to these tests. In the case of the deep model of HD 209458b, the predictions for the wind speeds are found to depend upon 
the chosen magnitude of the hyperviscosity, although the qualitative features of the simulated wind and temperature maps remain largely 
invariant. 

Within the spectral dynamical core, the (linear) fluid dynamical quantities are expressed as a truncated sum A'^h of spherical harmonics 
— the larger the value of A^h, the higher the horizontal resolution. The fiducial resolution we will adopt in the present study is T63 LA^v 
(A^h = 63), which corresponds to a horizontal grid of 192 by 96 points in longitude versus latitude. (See Table 2 of iHeng. Menou & PhillippsI 



^ This expression is considered to be order-of-magnitude because there is a factor of order unity associated with the exoplanetary albedo and the efficiency of 
heat redistribution from the day to the night side. 

^ The adjectives "shallow" and "deep" refer to the fact that one and about 5 orders of magnitudes in vertical pressure are simulated, respectively. 
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Table 1. Table of parameters, symbols and their values 



Symbol 


Description 


Units 


Earth 


Hot Jupiter 




stellar irradiation constant 


Wm-2 


938.4 


9.5 X 10^ 


At 


meridional temperature gradient parameter 


— 


1.4 




"s 


power law index for shortwave optical depth 




z. 


1 
1 


"L 


power law index for longwave optical depth 




4 


2 


•^So 


surface optical depth of shortwave absorbers 




0.2 


1401 




surface optical depth of longwave absorbers at equator 




6 


4.67 X 10^ 




surface optical depth of longwave absorbers at poles 




1.5 


4.67 X 10^ 


h 


strength of well-mixed longwave absorbers 




0.1 


1/2000 



cp specific heat capacity at constant pressure (of the atmosphere) J K ^ kg ^ 1004.64 14308.4 

Tl ideal gas constant (of the atmosphere) JK-^kg-^ 287.04 4593 

K, = H/cp adiabatic coefficient* JK-^kg-l 2/7 0.321 

Po reference pressure at bottom of simulation domain bar 1 220 

Cint areal heat capacity of lower atmospheric boundary JK-^m" 10^ 10^ 

acceleration due to gravity ms"^ 9.8 9.42 

Rp radius of (exo)planet km 6371 9.44 x 10* 
rate of rotation of (exo)planet s"! 7.292 X 10"^ 2.06 X lO'^ 





critical bulk Richardson number 




1 




^rough 


roughness length 


m 


3.21 X 10-5 




KvK 


von Karman constant 




0.4 




h 


surface layer fraction 




0.1 






hyperviscous time 


dayt 


0.1 


10-5 


At 


time step 


s 


1200 


120 


Tinit 


initial temperature 


K 


264 


1824 




vertical resolution 




20 


33 



<S> latitude degrees -90°-90° -90°-90° 

e longitude degrees 0—360° 0—360° 

To=gp/cp dry adiabatic lapse rate Kkm— ^ Ri 9.8 ~ 0.7 

6»T potential temperature K ss 300-500 ~ lO^-lO" 

'f Eulerian mean streamfunction kgs~^ 10^-10^" ~ lO^'^-lO^' 



f : expressed in terms of a (exo)planeta ry day (i.e., rotational period). 
|: we term k the "adiabatic coefficient", foUowing lPierrehumberll bOlOl) . and avoid calHng it 

the "adiabatic index" in order to not confuse it with the ratio of specific heat capacities. 
Note: one "Earth day" is exactly 86400 s (as used in the text and simulation/analysis codes). 



boil^ for a list of and their corresponding horizontal resolutions.) A finite difference solver with A'v grid points is used for the vertical 
coordinate, which is cast in terms of the "a-coordinate" ( iPhilliDSi.l957i) . 

(4) 

where P denotes the vertical pressure and Ps is the time-dependent pressure at the surface. By contrast, the reference pressure at the bottom 
of the simulation domain, Po, is a constant. Analogous to a, we define ctq = P/Po- For Earth-like simulations, the reference pressure is 
usually set to Po = 1 bar = 10^ Pa = lO'^ dyn cm"^. 

Numerical noise accumulates at the grid scale and has to b e damped via a "hy perviscous" term with a damping order of = 4, such 
that the operator V^"" acts on the relative vorticity (see §3.1 of iHeng. Menou & P hillipps 2011). The pragmatic aim is to apply horizontal 
dissipation on a time scale that is a small fraction of an exoplanetary day. For example, in the Held-Suarez benchmark for Earth, the 
a dopted time scale is ti, ~ 0. 1 Earth day. For our hot Jupiter simulations, we use ti, ~ 10-^ HD 209458b day, consistent with the value used 
in Heng. Menou & Phillipp sI (2011); this is also, to within an order of magnitude, the highest value needed. It is important to keep in mind 
that the application of hyperviscosity and the specification of is unsupported by any fundamental theory — horizontal dissipation is solely 
a numeric al tool meant to pre vent a simulation from failing due to spectral blocking. 

As in lHeng. Menou & Ph illipps ( 201 1), the simulations are started from an initial state of windless isothermality with Tinit denoting the 
initia l temperature, in contrast to using initial temperature-pressure profiles obtained from one-dimensional radiative-convective calculations 
(e.g.. lShowman et a l. '2009':'Dobbs-Dixon, G umming & Li jboiolJLewis et alj2O10l) . The possible dependence of our results on more exotic 
iiutisil conditions (Thrastarson & Cho 2010) is deferred to a future study. There is an initialization period for each simulation which we 
disregard, during which the temperatures and velocities, averaged over the entire (exo)planet, are decreasing/increasing from some initial 
values to their quasi-steady values. For Earth-like simulations, we find the initialization period to be 200 days. For hot Jupiters, a longer 
initialization period of 500 days is needed. After quasi-equilibrium is attained, we execute the simulations for a further 1000 days and 
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temporally average the ouput. Minor asymmetries in the physical quantities between the northern and southern hemispheres are artifacts of 
averaging over a finite period of time (i.e., in this case, 1000 days). 



2.2 Stellar irradiation and shortwave absorption 



At the top of the atmosphere, stellar irradiation is specified via the function 

-^TOA ~ ^oQ, 

where To is the stellar irradiation constantlf] 

' 1370 W m" 
9.5 X 10^ Wm 



a 



""SB-/* : 



R 



5780 K 



(1.146*^0) ( 0.0468 Au) (sooo k) 



(Earth), 

(HD 209458b), 



(5) 



(6) 



.146 Rq 

where ctsb denotes the Stefan-Boltzmann constant and the dimensionless stellar irradiation function Q — Q{^, O) is in general a function of 
both the latitude $ and the longitude Q. Since hot Jupiters orbit at distances o ~ 10 i?« from their host stars, we have To ~ 6.4 x lO"" W 
m^^ for i?* — Rq. By contrast, we have J-q ~ 51 W for Jupit er, which is a factor ~ 1 0"^ less than for hot Jupiters. In the case of HD 
209458b, we use the values of 7?^,, a and T.* measured or inferred bv lMazeh et al. I teoool) and lBrown et al.l (120011) to obtain the value shown 
in equation (|6f; we note that a/R^, ~ 9. In general, the influence of the stellar irradiation on a (exo)planet and its dependence on geometry 
(through Q) is known as the "thermal forcing' 



T he shortwave optical depth rs determines how the stellar irradiation is absorbed as it travels downward into the atmosphere (iFrierson 

2007al : lo'Gorman & Schneidej2008l:lMerlis & Schneidejboiol) , 



rs 



rSnf^O 



(7) 



Specif ying the power law index as ns 
201 Ibh 



is equivalent to having a uniformly mixed shortwave absorber, since jGuillofcoiol :! iHeng et al 



TS = = 10 

9p 



9p 



^0.01 cm2 g-i 100 bary V 10 m s-i ' 
where ks is the shortwave opacity and Qp is the surface gravity. When is constant, we have rs oc P. When ns > 1, the shortwave 
absorber is less well-mixed and resides lower down in the atmosphere. The received flux as a function of vertical pressure is then 

= Ttoa exp (-rs) . (9) 

Finally, we note that it is possible to specify a mean (exo)planetary albedo, but this only serves to diminish the value of J-toa. As such, 
one only needs to explore the effects of varying the stellar irradiation constant. 



2.3 Radiative transfer and longwave absorption 

The radiative transfer scheme is based on the two-stream approximation jMihalaJl978l : 



Hanse 



iTiooi) 



, which assumes that the radiation 



can be separated into longwave, upward- (Tf^) and downward-propagating (^^.^l) fluxes (see §4.2 of .Pierrehumbert.2010i) . 



drL 



T 



tL 



(10) 



+ J"b 



(bottom) 



where T^h = os-rT is the blackbody flux. At the bottom of the simulated atmosphere, the longwave boundary condition is T^_^ — 
o'ssTg where Ts is the surface tempe rature; at the top, it is -^[|^° ''^ = 0. The Newtonian relaxation term in the temperature equation is 
replaced by a radiative source term (iFrierson. Held & Zurita-GotoilboOdI) , 

Qr.d = ~^-^{J't^-J'i^-J'is), (11) 

where cp is the specific heat capacity at constant pressure, p is the mass density and z is the vertical coordinate. The shortwave, upward- 
propagating flux is assumed to be zero, i.e., = 0. Furthermore, substituting the Newtonian relaxati on scheme for the dual-band ra- 
diati ve transfer scheme alleviates t he need to rely on one-dimensional radiative-convective models (e.g., Ilro. Bezard & Guillotll2005l : see 
also Ramanathan & Coaklev 19781) for calculations of the equilibrium temperature-pressure profile, thereby allowing self-consistency to 
be achieved. From the perspective of equation (110) alone, the radiative transfer scheme is strictly speaking grey (as already noted by 
Frierson. Held & Zurita-Gotorl2006 ^ — the shortwave, downward-propagating flux is merely diluted on its way down into the atmosphere via 



equation ^ — but we term it "dual-band" to emphasize the dual-wavelength nature of our treatment. The effects of scattering are neglected. 



^ The Solar constant is typically taken to be about 1367 W m ^ 
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The bottom of the atmosphere (or the "surface") is an ideahzed slab with a specifiable heat capacity, described by a single temperature 
Ts governed by the equation, 



(12) 



where Cint is the (areal) heat capacity of the slab. Horizontal transport of fluid within the slab is not modelled. In the (terrestrial) atmospheric 
sciences literature, this idealized slab mimics the "mixed layer ocean", which is the ~ 100 m-thick layer of ocean inteimediate between 
the atmosphere and the deep ocean, characterized by rigorous mixing induced by wave motion and turbulence. For simplicity, we do not 
implement an internal heat flux, typically parametrized by the equivalent blackbody tem perature Tint. For example, Ir o. Bezard & GuiUotl 
JioOS) and Guillot (2010) use Tint = 100, 300 K in their models of HD 209458b, while Liu & Schneid er (2010) adopt 5.7 W m'^ as the 
internal heat flux in their simulation of Jupiter (which is equivalent to Tint ~ 100 K). 

As the radiation is absorbed and re-emitted, the (infrared) photons encounter a longwave optical depth near the surface, described by 

(13) 



-TLn 



where the longwave optical depths at the equator and the poles are given by tl^^ and TLp^,^^ , respectively. As already noted bv lMerlis & Schneider 
(20l3), equation l ll3t ignores the longitudinal dependence of introduced by water vapour feedback and is thus inadequate for treating 
tidally locked. Earth-like aquaplanets. Throughout the vertical extent of the atmosphere, the longwave optical depth is a linear combination 
of well-mixed and segregated atmospheric absorbers, 

tl = tl„ (To + , (14) 

where rL„ and tl^ generally depend on geometry (i.e., $ and Q) as well as atmospheric chemistiy The second term in equation ( 114b 
represents species of absorbers where the associated pressure scale height is I/wl that of the well-mixed species. For example, the lin- 
ear term may represen t well-mixed species like carbon dio xide, while setting wl = 4 approximates the structure of water vapour in the 
terrestrial atmosphere ( Frierson. Held & Zurita-Gotor 20061) . As another example, setting riL = 2 app roximates the process of collision- 
induced absorp tion ('Herzberg^'l952') within the atmospheres of the giant planets in our Solar System ( Liu & Schneider 2010h . Following 
Frierson. Held & Zurita-Gotor, (.2006> .2007.) and Frierson ( 2007 ah . we use 

TLw = T"Lo/i, 
TL, = TLo (1 - }l) , 

with fi being a dimensionless parameter that contro l s the strength of the well-mixed longwave absorber. 



(15) 



In 



Frierson. Held & Zurita-Gotod (l2006ll2007n . 



Frierson 



(l2007ij) and lO'Gorman & Schneiden ( 1200 8h . a key assumption is that the (v 



ter) moisture content does not affect radiative transfer, such that th e dynamical effects of the latent heat release associated with the con- 
densation of water vapour can be isolated from its radiative effects. iMerlis & Schneider ( I2OI0I) improved upon this treatment by allowing 
the optical depth associated with water vapour (tlj, ; wl =4) to depend on the instantaneous column density, thereby providing a simple 
treatment of water vapour feedback. Since we are concerned only with dry atmospheres in the present study, we will ignore this issue. 



2.4 Convection 



In simulations of atm ospheric circulation, convection is typica lly mimi cked using simplified parametri zations known as "convective 
adjustment" schemes (e.g., Manabe, Smagorinskv & Stricklej|l965 : §3.6.9 of Washington & Parkinson 2005). For a dry atmo sphere in hy- 
drostatic equilibrium, the change in temperature of a parcel of atmosphere resulting from vertical motion is (e.g., §2.7.2 of ltlolton 2 004^1 



where the lapse rate is defined as 



1 OBt 1 



dT 



(16) 



(17) 



with z being the vertical spatial coordinate, and the dry adiabatic lapse rate is To = Qp/cp, which has a value of about 9.8 K km ^ for the 
terrestrial atmosphere. The potential temperature, 

0T = Tao^, (18) 

where k = TZ/cp and TZ is the ideal gas constant, is the temperature a parcel of atmosphere would have if it was compressed or expanded 
adiabatically from the pressure P to the reference pressure Pq at the bottom of the simulation d omain. The adiabatic coefficient k is related to 
the number of degrees of freedom of the atmospheric gas ndof : k ~ 2/(ndof + 2) (see §2.3.3 of |Pierrehumbertl2010t) . For example, molecular 
hydrogen is experimentally found to have ndof ~ 5.2, close to the theoretical value for a diatomic gas (ridof = 5). The atmosphere is unstable 
if F > To (i.e., super-adiabatic lapse rates)lf| When the lapse rate is adjusted to a state of stability — subjected to the constraint of total energy 
conservation — convective adjustment is said to be performed. In other words, the tendency of convection is to minimize the magnitude of 



* Equation (T6j may be generalized to describe an atmosphere containing moisture and with a corresponding moist adiabatic lapse rate. 
^ The condition for stability is sometimes known as the Schwarzschild criterion. 
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the po tential temperature gradient OOt /dz and convective adjustment is an approximate way to mimic this process jRamanathan & Coaklev 



1978 



Within the FMS, the dry convective adjustment scheme follows the prescription of lManabe. Smagorinskv & Stricklen(ll965h . At a given 
time, the temperature difference 5T resulting from departures from the adiabatic lapse rate is computed, 

dP 



= 0, 



— I STdP - 
9p 



(19) 



0, 



where the integral is taken over the vertical extent of the unstable, model atmosphere layer. An implicit assumption made is that when the 
lapse rate of a layer is super-adiabatic, convection is rigorous enough to maintain a neutral lapse rate of potential temperature. The kinetic 
energy created by convection is then dissipated and instantly converted into heat, such that the total potential energy is invariant to convective 
adjustment. Within each model layer, the temperature is then adjusted by ST. 

It is worth noting that there exists a body of work on moist convective adjust ment schemes, largely motivated by the pioneering work of 



Manabe. Smagorinskv & Strickleii ( ll965n . |Betts'('l986') and'Setts & Millejjl986h . who proposed the convective adjustment to reference tem- 



perature and humidity profiles instead of Fq . Lprierson ( 2007a,) adopted a simplified version of the Betts-Miller scheme, where the reference 
humidity profile is held at a constant threshold value. 



2.5 Atmospheric boundary layer 



The boundary layer is the part of the atmosphere where the flow is strongly influenced by interaction with the surface. It can be thought 
of as enforcing a "no slip" boundary condition, s uch that large velocity gradients — and hence shears — exist across a small vertical 
height, which induces turbulent mixing (e.g., §5 of lHoltorJ l |2004) . The boundary layer may become convectively unstable, which induces 
mixing. The end result is the same: eddies are produced which are smaller in size than those resulting from large-scale rotational flows, 
are typically unresolv ed in global simulations of atmospheric circulation and thus need to be accounted for using sub-grid models (e.g., 
Troen&Mahr3ll986h . In other words turbulence within the boundary layer may be driven mechanically or via buoyancy. The eddies tend 
to have comparable horizontal and vertical length scales — i.e., they are effectively three-dimensional — and are responsible for heat and 
momentum exchanges between the atmosphere and the surface. Consequently, the boundary layer affects the d ynamics and th ermodynamics 
of the terrestrial atmosphere and is responsible for a non-negligible fraction of kinetic energy loss by the flow ( Garrattiri994 ). On Earth, the 
height of the boundary layer is ~ 30-3000 m and contains ~ 10% of the mass of the terrestrial atmosphere. 

We find that the boundaiy layer scheme is needed if the bottom of the simulation domain is placed within the active part of the 
atmosphere, such as in an Earth-like simulation. However, this scheme is not required to successfully execute hot Jupiter simulations. 
Moreover, it is not entirely clear if the analogy with a terrestrial boundary layer can be carried over to mimic drag in the deep, inert layers 
of a hot Jupiter. As such, although we include a description of our boundary layer scheme for completeness (since it is used for our Earth- 
like simulations), we relegate the technical details to Appendix |A] As described therein, our atmospheric boundary layer scheme requires 
the specification of four additional parameters: the roughness length Zrough, the critical bulk Richardson number 7?.,;,crit, the von Karman 
constant KvK and the surface layer fraction /;,. 



2.6 Additional physical quantities 



In the interest of clarity, we will define physical quantities, used in the present study, which may be unfamiliar within the astronomi- 

cal/astrophysical literature. 

The potential temperature 6t (equation 1181 ) is related to the entropy 5 by (e.g., §1.6.1 of Vallisl20o3 . equation [3.7] of Peixoto & Port 



1984) 



5 = cp In 6t- 



(20) 



The preceding expression is only valid if cp is constant (otherwise it involves an integration over 6t)- Therefore, examining profiles of the 
potential temperature is equivalent to determining where the surfaces of constant entropy ("isentropes") reside, which provides insight into 
the convective stability of the atmosphere. 

Insig ht into the large-scale circulation patterns within the atmosphere may be obtained by examining the Euleran mean streamfunction, 
defined as jPeixoto & o'or]|l984l:|Pauluis. Czaia & Kortvll200^ 



RpPo cos<l? 



(<&, CTo) dao 



(21) 



where vg, — v,s>{$, ao) is the temporally- and zonally-averaged meridional velocity. On Earth, we have ~ 10® kg s^^. We note that the 
"Sverdmp" (Sv), where 1 Sv = 10^ kg s~^, is an alternative unit used in t he oceanography community. I n the case of hot Jupiters we 
have 'if ~ 10^* kg s~^. Our definition of the streamfunction follows that of lPauluis, Czaia & Kortvl (20081) , such that ^' < for a circulation 
cell situated just above the equator in the northern hemisphere, there is southward flow at low altitudes (high P values) and northward flow 



© 2013 RAS, MNRAS 0Q0.[Tll28] 



Atmospheric circulation of exoplanets 11 7 



Held-Suorez benchmark 



Held — Suorez benchmark 





Figure 1. Temporally-averaged, zonal-mean zonal wind (left column) and temperature (right column) profiles as functions of vertical pressure P and latitude 
Top row: Held-Suarez dynamical benchmark. Middle row: Earth-like model without convection. Bottom row: Earth-like model with convection. Contours 
are in units of m s~^ (left column) and K (right column). 



at high altitudes (low P values). It is worth noting that the presence of vertical motions due to atmospheric circulation does not violate 
the assumption of hy drostatic balance provided the vertical accelerations are small compared to the acceleration due to gravity (see §2.5 of 
Pierrehumber3l2O10l) . 

In Figures l6l [8l 1 1 21 and 1 141 we opt to split up the Eulerian mean streamfunction into the day and night sides of the exoplanet. Strictly 
speaking, this implies that our use of the word "streamfunction" to describe the relevant panels in these figures becomes invalid, since 
streamfunctions are normally associated with circulations which are non-divergent in the plane being considered — an integration over all 
longitudes is strictly required. However, our intention is to capture the differences between the atmospheric circulation on the day and night 
sides of a tidally locked exoplanet. The reader should therefore be mindful of the way we use the term "streamfunction" in these figures. 
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Held-Suorez benchmark 



Held — Suarez benchmark 





Figure 2. Temporally-averaged, zonal-mean potential temperature (left column) and the Eulerian mean streamfunction (right column) profiles as functions of 
vertical pressure P and latitude <E>. Top row: Held-Suai'ez dynamical benchmark. Middle row: Earth-like model without convection. Bottom row: Earth-like 
model with convection. Contours are in units of K (left column) and X lO" kg s~^ (right column). 

3 EARTH-LIKE MODELS 



3.1 Earth 

The first step is to esta blish a Earth-like model, both as a consistency check and as an operational baseline from which to generalize 
to tidally-locked e xoplanets jHeng. Menou & PhilUppi [201 1 : Heng & Vogt 2011 ). The present model should be regarded as a dry, simpli- 
fied version of the Iprierson. Held & Zurita-Gotor ( 20061) model. We omit latent heating effects fo r simplicity and thus the (dry) Earth-like 
models presented here have less fidelity in simulating the terrestrial climate than the models of Frierson. Held & Zurita-Goto3 (2006) in 
aspects which are strongly influence d by water condensation, such as lapse rates and meridional overturning circulation strengths. The 
Frierson. Held & Zurita-Goton (120061) model with moisture has a significantly better fit to observations than the Held-Suarez model in as- 
pects of the terrestrial circulation such as the strength of the Ferrel cell, meridional energy transports and radiative cooling rates. Following 
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Figure 3. Globally-averaged temperature-pressure profiles for the trio of Earth-like simulations presented in ^3.1\ 



Frierson. Held & Zurita-Gotoi ( 20061 2007 ) and Frierson ( 2007al) . we define the stellar irradiation function to be 

At 



5 = i 



1 + — ^(l-3sin^ $) 



(22) 



ine ioiar constant is tauen to oe j-q ~ yo8.4 w m , wnicn corresponas to a giooai-mean aioeao oi aoout jU7c. ine seasonal cy 
insolatior[f| is not considered, meaning the exoplanet is implicitly assumed to reside on an orbit with zero eccentrici ty and obliquity , 
aside, we note that solar irradiation has been inferred to vary ove r much longer tim e scales (Berger & Loutre 1991; Crowley 2000j. 



The meridional temperature gradient can be adjusted by varying the value of the dimensionless constant At, which we take to be 1.4. 
The Solarconstant is taken to be J"o ~ 938.4 W m^^, which corresponds to a global-mean albedo of about 30%. The seasonal cycle of 

As an 

longer tim e scales (Berger & Loutre"199r; 'Crowlev'2000^. Other 
parameter values are listed in Table[T]and are largely adopted from lFrierson. Held & Zurita-Gotor (2006) and Frierson. (2007a ). 

A key difference between the Held-Suarez model and our dry Earth-like model is that the radiative equilibrium of the latter is strongly 
unstable to convection. Preliminary simulations without the boundary layer scheme ([|23} implemented fail, suggesting that models in which 
the surface is located within the active — as opposed to inert — part of the atmosphere require a scheme to establish heat and momentum 
equilibrium between the atmosphere and the surface. This expectation is borne out in our hot Jupiter simulations — where the surface is 
located at P ~ 100 bar, well into the inert part of the atmosphere — which do not require the boundary layer scheme to run to completion. 
We therefore include the boundary layer schemes only in our Earth-like models. 

In Figure [T] we show the temporally- and zonally-averaged profiles of zo nal wind and tempera ture as functions of the vertical pressure 
P and the latitude $, which we term the "Held-Suarez mean flow quantities" (Held & Suarez|[l994l ). For comparison, we include the Held- 
Suarez dynamical benchmark, as computed by Heng. Men ou & Phi llipps (2011). In the trio of simulations, the jet streams in the upper 
troposphere (P ~ 0.3-0.4 bar), at mid-latitudes, are clearly evident. There are quantitative differences between each of the temperature 
profiles, but the general trends are that the temperature decreases away from the equator and also with increasing altitude. 

More insight on the convective stability of the model atmospheres may be gleaned by examining the temporally- and zonally-averaged 
po tential temperature pro files, as shown in Figure |2] The potential temperature may be regarded as the "height-adjusted" temperature (§111 
of |Peix6to&Oortlll984l) . If it monotonically increases with height and does not depend on latitude, its structure is termed "barotropic". 
In this case, no convection occurs. If there is a gradient in the potential temperature across latitude, then its structure is said to be "baro- 
clinic". In suc h situations, horizontal or "slanted" convective motions occur Tapping into the available potential energy of the atmosphere 
Lorenz'1955'), he baroclinic eddies resulting from these motions are responsible for transporting sensible heat from the equator to the poles 
Pier rehumbert & Swanson 1995). The stratosphere and troposphere are thus regions of the atmosphere where the potential temperature 
structure is (largely) barotropic and baroclinic, respectively, at least within the temporally- and zonally-averaged context of our models. In 
reality, the terrestrial stratosphere is not strictly barotropic — along isobars, there are distinctly non-zero variations of temperature with 
latitude which are dynamically important. The h eight or vertical pressure at which the troposphere transitions into the stratosphere is the 
tropopause (e.g., see Figure 4 of Frierson|[2007a ), located at P ~ 0.2 bar near the equator (the tropics) and P ~ 0.3-0.4 bar otherwise (the 
temperate and polar regions). 

The middle row of Figure |2]shows the simulation with radiative transfer but without convective adjustment. Unlike in the Held-Suarez 
benchmark, the existence of a surface with a finite heat capacity (see i]2.3b and its heating by solar irradiation (see i]2.2b results in con- 
vectively unstable vertical columns, thus necessitating the treatment of convection. The bottom row of Figure [^demonstrates how the dry 
convective adjustment scheme (see i]2.4b brings the model atmosphere to convective stability in the vertical direction, while preserving the 
baroclinic instabilities (i.e., horizontal convection) produced in the simulation. Essentially, the role of convective adjustment is to straighten 
the isentropes such that they are vertically stable. 



In this study, we will use the terms "irradiation" and "insolation" interchangably. 
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Tidally-locked Earth (without convection) 



Tidolly-locked Earth 




Tidally-locked Earth (without convection) 




Tidally-locked Earth (without convection) 





Figure 4. Held-Suai'ez mean flow quantities for the tidally-locked Earth simulation. Top row: zonal wind (m s ^). Middle row: temperature (K). Bottom row: 
potential temperature (K). Left column: without convection. Right column: with convection. 



Atmosphe ric circulation of the terrestrial a tmosphere can be generally regarded as having a three-cell structure (e.g.. |Peix6to & Port 

1984 §2 .1.5 of IWashington & Parkinsonll2005h . The Hadley and polar cells, located near the equator and poles, respectively, are direct 
circulation cells because they are driven by heating and cooling patterns of the Earth's surface. Both cells are characterized by rising air in 
their warmer branches. By contrast, the mid-latitudinal Ferrel cells are indirect because they are driven by the presence of baroclinic ed dies; 
cooler air is forced to rise. The p resence of these cells may be revealed by examining the Eulerian mean streamfunction (e.g., Frie rsonI 
2007al : iMeriis & Schneide^boid) . as defined in equation ( I2U . The Held-Suarez benchmark, which mimics heating and cooling patterns 
by forcing the temperature field to relax to an ad hoc "equilibrium temperature" profile on a specified "cooling" time scale, manages to 
produce the Hadley, Ferrel and polar circulation cells, although the circulation associated with the Ferrel cells is too strong by a factor of two 
(top right panel of Figure |2). As observed on Earth, the polar cells are relatively weak with *!/ having values about an order of magnitude 
lower than for the other two cells. The simulations with dual-band radiation (middle right and bottom right panels of Figure |2j manage 
to only produce the Hadley and Ferrel cells; there is a pair of cells near the poles, but they reside only at low altitudes {P ~ 0.9-1 bar). 
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Figure 5. Wind maps from the tidally-locked Earth simulation with convective adjustment. Left column: zonal wind. Right column: meridional wind. The 
first, second and third rows show the maps at P = 0.25, 0.55 and 1 bar, respectively. Contours are in units of m s~^. 



The absence of the polar cells in our simulations and also those of iFrierson. Held & Zurita-Gotor (I2OO6I) is unsurprising since the details 
concerning the polar temperatures (e.g., sea ice, the presence of Antarctica) are not modelled. The Hadley cell is significantly too strong in 
the dual-band simulation due to the near-neutral convective stability in the tro pics, while the Ferrel cell is confined to only high altitudes due 
to the presence of convective momentum fluxes in the lower troposphere (see lprierson. Held & Zurita-Gotorll2006 for more details) — these 
issues are ameliorated in simulations with moisture. 

Figure[3]shows the globally-averaged temperature-pressure profiles for the trio of simulations. By construction, the Held-Suarez bench- 
mark does not produce the temperature inversion empirically observed in the stratosphere, since the stratospheric temperature is set to be 
constant in this model. There are slight differences between the T-P profiles for the simulations with and without convective adjustment, 
where the latter model produces a small temperature inversion in the stratosphere. 

Our dry model neglects the fact that the terrestrial tropospheric temperature profile is strongly in fluenced by latent hea t release from 
the condensation of water vapour. In the tropics, the temperatures are approximately moist adiabatic (IXu & Emanuell[l989i) . while in the 
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Tidally-locked Earth (without convection) (day side) 



Tidolly-locked Earth (day side) 





Figure 6. Held-Suarez mean flow quantities for the Euleiian mean streamfunction (in units of xlO'"* kg s^^) in tlie case of the hypothetical tidally-locked 
Earth. Top row: day side. Bottom row: night side. Left column: without convection. Right column: with convection. The irregular contour inter'vals are chosen 
to reveal the presence of weaker circulation cells. 



extratropics they are determined by a coupling between latent heat release and baroclinic eddies (e.g., Frierson. Held & Zurita-Gotoijbood) 
with additional contributions from the land-sea contrast. On average, these factors assure that the tropospheric lapse rate on Earth is less than 
the dry adiabatic value. 

In summary, we have generalized our Earth-like model away from the purely dynamical Held-Suarez benchmark and demonstrated the 
utility of examining profiles of the potential temperature and Eulerian mean streamfunction. We have also demonstrated the necessity of 
including the treatment of convection alongside radiative transfer, at least for our Earth-like models. We next extend the same simulation and 
analysis techniques to a hypothetical tidally-locked Earth. 



3.2 Tidally-locked Earth 



The case of a hypothetical tidally-locked E arth is a useful and computation ally efficient case study for operat ionally transitioning 
between the simulation of Earth and hot Jupiters jHeng. Menou & PhillippsI 201 1 ). The stellar irradiation function is ( Merlis & Schneider 
20ld) . 



e = max{0,cos$cos(e-eo)}, (23) 

where the substellar point is located at $ = 0° and Q = Qq. Our (arbitrary) choice is Qq — 180°. The parameter values are the same as 
thos e adopted in the Earth-like simulation, but the rotational freq uency is slowed down such that the exoplanet takes one year to rotate on its 
axis jMerUs & Schneidejboid iHeng. Menou & Phillippj201lh . 

Qp np/365. (24) 

Our functional form for the normalizati on of the longwave optical depth, as descr ibed by equation il3\, does no t have a longitudinal 
depe ndence and is the same as that adopted bv lFrierson. Held & Zurita-G otor ( 2006) and O' Gorman & Schneider ( 2008 ). Merlis & Schneidei 
mid) have noted that a more realistic treatment is to allow for tlq to longitudinally vary due to longwave water vapour feedback. However, 
since our intention is to implement a simplistic, dry model for a tidally locked exo-Earth merely as a sanity check, we ignore this refinement 
while acknowledging its need if the intention is to model tidally locked. Earth-like aquaplanets. 
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The Held-Suarez mean flow quantities for zonal wind, temperature and pote ntial temperature are shown in Figure 2] Our results for 
the zonal-mean zonal wind profile are similar to the top-left panel of Figure 6 of iMerlis & Schneidej fcoid), who noted that the profile 
contains weak residuals of opposing contributions from various longitudes. It is therefore unsurprisingly that the simulations with and without 
convective adjustment have zonal-mean zonal wind profiles which are somewhat different. This phenomenon can be seen more clearly in the 
zonal and meridional wind maps as functions of longitude and latitude (Figure|5j, which indicate the presence of large circulation cells. The 
potential temperature profile shows that the atmosphere is, on average, approximately barotropic down to P ~ 0.5 bar, implying that the 
tropopause is located farther down. We show the zonal-mean potential temperatures without separating them into day and night side profiles 
as they are fairly similar; we will examine this issue in greater detail for the hot Jupiter simulations ([Q. Rather, our main intention is to 
demonstrate the effect of convective adjustment. 

The Eulerian mean streamfunction in Figure |6]reveals the presence of multiple circulation cells on the day and night sides. On the day 
side, a pair of cells exist in both the northern and southern hemispheres. The cells are able to extend from the equator to the poles, because 
the reduced rate of rotation of the exoplanet weakens the Coriolis deflection. On Earth, the relatively rapid rotation limits the Hadley cells 
to about c[> ^ ±30°. On the night side, a pair of large, direct cells reside at higher altitudes (P < 0.6 bar); a collection of four smaller 
circulation cells is evident at lower altitudes. These features are robust to the application of convective adjustment. 

The hemispherically-averaged temperature-pressure profiles for the day and night sides of the tidally-locked Earth are shown in Figure 
[T] On average, the temperature differences between the two sides are barely 10 K and are only evident at P > 0.6 bar As in the Earth-like 
case (i )3.1b . convective adjustment has little effect on the T-P profiles. 



3.3 Other studies 



The computatio nal setup descri bed here has also been used for other studies of the terr estrial atmosphere, includin g the analysis of 
equatorial transients ( Friersorj2007bb. the width of the Hadley circulation in different climates ( Frierson. Lu & Cherj2007 ), the vertical tem- 
perature profile in mid-latitudes (^Frierson'200 8|), how the cross-equatorial Hadley circulation is affected by remote heating /Kang et al.'2 008h 
and t he location of the jet streams (iLu. Chen & Friersonli201ft) . It was also applied to the study of the atmosphere on Titan (iMitchell et alJ 
2009h . 



4 HOT JUPITER 



To facilitate comparison with our prev ious work, most of the parameter valu es we adopt for our hot Jupiter sim ulati ons are culled from 

the de ep model of HD 209458b presented in lCooper & ShowmanI J2005I.1 20061) and lHeng. Menou & PhilUppsI hoilh . As in lHeng. Menou & Phillipps 
1 20111) , we use a T63L33 resolution with a range of pressures of 1 mbar < -P ^ 220 bar. A noteworthy difference from the Earth-like models 
is the assumption that stellar irradiation impinges upon a well-mixed shortwave absorber of unspecified chemistry (ns = 1). The adiabatic 
coefficient used is k — 0.321 (ridof ~ 4.2). 



4.1 Previous work: purely dynamical model 

As a prelude to our resu l ts, we re-analyze the simulation output from the purely dynamical, deep model of HD 209458b as computed 
by Heng. Menou & Phillipp sI i 2011 ). Figure [S] shows the Held-Suarez mean flow quantities for the potential temperature (using equation 
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Figure 8. Held-Suarez mean flow quantities for the T63L33 deep model of HD 209458b (see text). Top row: potential temperature (K). Bottom row: Eulerian 
mean streamfunction (x 10^^ kg s^^). Left column: day side. Right column: night side. For the potential temperature profile, contour intervals are made uneven 
in the interest of clarity. For the Eulerian mean streamfunction profile, one contour interval is irregular so as to reveal the presence of weaker streamlines. 



nisi ) and Eulerian mean streamfunction (using equation fllV ) separated into the day and night sides of the hot Jupiter. Near the top of the 
atmosphere (P ~ 1-10 mbar) on the day side, the atmosphere is slightly baroclinic but with flow converging towards the equator (unlike 
in the case of the terrestrial troposphere). On the night side, the atmosphere is noticeably more baroclinic. On both the day and night sides, 
the atmosphere becomes barotropic for P ^ 10 bar by construction, since the Newtonian relaxation time is defined to be infinite at these 
pressures. The relative weakness of Coriolis deflection again aflows for a pair of circulation cefl extending from the equator to the poles, as 
indicated by the Eulerian mean streamfunction. Unlike in the Earth-like (i j3.lt and tidally-locked Earth ( ij3.2t simula tions, the atmosphere 
sinks at the e quator and rises at the poles on the day side; the reverse happens on the night side. As demonstrated by Showman & Polvaiil 
1201 11) using a hierarchy of analytical models and simulations, there is downward transport of eddy momentum at the equators of 
hot Jovian atmospheres, which causes counter-rotating, equatorial flow. However, the ho rizontal convergence of ed dy momentum yields 
super-rotating, equatorial flow and tends to dominate the net flow (as shown in Figure 4 of Showman & PolvMiilboil ). Our simulations are 
consistent with this picture. It is somewhat counter-intuitive that the circulation is stronger on the night side (i.e., higher values of ^). There 
is virtually no meridional circulation in the inert {P ^ 10 bar) part of the atmosphere. 



4.2 Setup 



A number of con ceptual challenges exist when generalizing our setup to simulate hot Jovian atmospheres. Firstly, as discussed in 
Showman et al. I hoogh . it is presently unclear how to physically describe the frictional dissipation resulting from the differences in flow 
structures between the atmosphere and its deeper counterparts within a hot Jupiter. Therefore, while the boundary layer scheme described 
in i ]2.5l offers a tempting option to simulate frictional dissipation, we choose to switch it off for our hot Jupiter simu lations. Essentially, 
our models assumes "free slip" lower boundary conditions. A different approach is adopted by Liu & Schneider ( 2O10h . who simulated the 
zonal wind structures of Jupiter, Saturn, Uranus and Neptune by varying the (Rayleigh) drag at the lower boundary in order to better fit 
observations, with the justification that such a formulation mimics the Ohmic dissipation caused by the induced magnetic fields opposing the 
zonal flows. Since there are virtually no empir ical constraints on the z onal wind profiles of hot Jupiters, unlike for the giant planets in our 
Solar System, we do not adopt the approach of Liu & Schneidei ( 201(]|) . 
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Secondly, the bottom of the simulation domain ( q2.3b does not lend itself to easy physical interpre tation, unlike in the Earth- like cases 
where it mimics the mixed layer ocean. A plausible approach is to demand that the radiative time scale ( Showman & Guillo3l2002 ). 

-3 



^rad 



CpP 



350 days 



9.42 m s- 



r 



1700 K 



(25) 



^14308.4 J kg-^ 220 bai- 
ls continuous at the bottom. At and beneath the bottom of the simulation domain, the hot Jupiter has a finite thermal inertia if the heat capacity 
is non-zero (Cint 7^ J K 

2 



^ m ^) — both longwave and intrinsic heat are not instantaneously (re-)emitted. The thermal inertia time scale is 



tx 



Cint 



(26) 



where the thermal inertia is defined as (e.g., |Palluconi & Kieffeilll98lb 

X = (fccon Pint Cp^int)^^^ , (27) 

with fccon, Pint and cp,int being the themal conductivity, mass density and specific heat capacity at constant pressure associated with the 
bottom, respectivel y. The thermal conductivity ranges from ~ 0.1 WK~-^m'^ for gaseous hyd rogen (and helium) to ~ lO'^ W K ^ m ^ for 
metal lic hydrogen l lHubbardll968 ): the latter is not expected to exist below pressures ~ lO** bar JstevensonI 19821 ; iHubbard, Burrows & Luninel 



2002). 



To ensure continuity between the simulated atmosphere and the bottom, we dem and that cp,int 
eliminated in favour of P and T via the ideal gas law (e.g., §2.3.1 of Pierrehumbertlbo 1 ol) . 
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3 X 10^ gcm~^ 
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4593 J K"^ kg-i 1700 K 



Cp. Furthermore, pint may be 



(28) 



V220 bar^ 

where the specific gas constant isTZ = TZ* / fi, TZ* — 8314.5 J K^^ kg^^ is the universal gas constant and fi is the mean molecular weight 
of the atmospheric molecules. The specific gas constant may be related to more familiar quantities via 

k-R 



(29) 



where fee is the Boltzmann constant, m — /imn is the mean mass of the atmospheric molecules and niu denotes the mass of a hydrogen 
atom. For example, if 7?. = 4593 J K^^ kg^^ (as we are assuming for HD 209458b), then /i ~ 1.8, close to the value for an atmosphere 
dominated by molecular hydrogen. By equating equations ( 125 l l and l l26t , we obtain a plausible estimate for the thermal inertia, 
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Besides the order-of-magnitude nature of the estimate, the main uncertainty lies in the value to assume for fc, 
results are insensitive to these uncertainties. 



(30) 

con. We will see later that our 



4.3 Analytical formalism for temperature-pressure profiles: generalization of Guillot (2010) 



The most effici ent way to initiate a given simu l ation is with a tempera tur e-pressure profi le in radiative equilibrium, as may be computed 
using the models of Hubenv. Burrows & Sudarskv ( 2003 ). Hansenl i 20081) or lOuillotlhoid) . We simplify our initial condition even further: 
we assu me a constant tem perature Ti nit which does n ot depend on latitude, longitude or pressure. Since the radiative time scale increases with 
pressure ( IShowman & Guillot 20021 : llro. Bezard & Guillot 2005). the temperatures at low pressures rapidly equilibrate to values consistent 
with dynamical and radiative equilibrium, whereas the temperatures at depth (Too) re main close to ra diative equilibrium. It is therefore 
plausible to select Tinit = Too. In this sub-section, we generalize the analytical results of Guillot 1 201oh to obtain the temperature-pressure 
profile with both til = 1 and 2, and subsequently use it to estimate the values for Too- 

By retaining th e assumption o f a constant shortwave opacity ks, but allowing for an arbitrary functional form for the longwave opacity 
kl, equation (29) of lGuillol may be generalized as 



3^4 
^Jint 



+ / KL dx' + 3/r, 



exp (^-~\/3ksx^ + ^ kl exp ^— \/3ks2;'^ 



dx' 



(31) 



where Tint is the blackbody temperature of the internal heat flux, kl ~ hii^{x'), dx = pdz, p is the mass density, z — Rp — z and z is the 
vertical spatial variable (with the centre of the exoplanet as a zero reference point). The distance from z = to the top of the atmosphere is 
Rp; we expect Rp ~ _Rp. The quantity x is thus the column mass, per unit area, measured from the top of the atmosphere. 

The quantity / ~ 0.25 is a dilution factor meant to mimic stellar irradiation being the strongest at the substellar point, since equation 
( I3II > assumes isotropic irradiation. The equilibrium temperature of the hot Jupiter, assuming zero albedo and no heat redistribution, is 

1/2 / T \ 1/4 



Tea = ^ 



T* 



V4crsB J 



(32) 
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Figure 9. Analytical temperature-pressure profiles calculated using our formalism in ^43] (solid curves), compared against simulated ones with only one 
Earth day of integration (dotted curves). Note that the simulated profiles are hemispherically-averaged over the day side only. Left panel: the increases in the 
temperatures at depth, due to the effect of collision-induced absorption, are described by a single parameter e. Right panel: varying the ratio of shoilwave to 
longwave opacity normahzations 70 for a fixed value of e = 2000. 



For the values of 7?*, a and associated with HD 209458b (see equation (SJ), we have Tcq ~ 1432 K. The corresponding irradiation 
temperature is Tirr = v^Tcq ~ 2025 K. 



If kl is constant, then equation | |3U reduces to equation (29) of Guilloll ( 2O10h . With Tint = K, it follows that the temperature at depth 



where 



2129 K/ 



1/4 



(33) 



(34) 



Note that the equality in the preceding definition only holds for the assumption of constant opacities. Equation (49) of Guillol 1 2O10l) . which 
averages over latitude and longitude, provides a more accurate estimate for the (mean) temperature at depthQ 



1539 K. 



Equating ([33} and J35I ) all ows for an estimate of the dilution factor required: / ~ 0.273. Following Guillot 1 2O10l) . we have adopted 7 — 0.6. 

As noted by .Guillol 1 2O1 0l). the assumption of til = 1 becomes inadequate deep inside the atmosphere, because collision-induced 
absorption becomes a non-negligible effect. The longwave opacity then scales linearly as pressure, which implies that a correction term with 
riL = 2 needs to be included. We demand that 

TLw + = TLo = ero when P = Po ox ao = 1, 

(36) 

TL„ ~ TQ when P <C Po or (To < 1, 

whence 

^' = ^' (37) 

TLg = (e — 1) To. 

The quantity to is the normalization for the longwave optical depth in the absence of collision-induced absorption, while the dimensionless, 
multiplicative factor e ^ 1 accounts for the increase of the longwave optical depth at P = Po due to this effect. Collision-induced absorption 
becomes important when 

Po 



We write the longwave opacity as 



P > 



«L = KO [1 -f 2 (e - 1) (To] 



KO 



l + 2(e - 1) 



where ko = QpTo/Po, xo = Po/Op and the second equality follows from the assumption of hydrostatic equilibrium, 

dP 



-pr = P9p 

az 



P = XQp. 



(38) 



(39) 



(40) 



^ If / = 0.25, equation (33) yields T^a « 1506 K. 
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Figure 10. Comparing analytical and simulated temperature-pressure profiles for the ng = riL = 1 hot Jupiter models. The profiles labelled "global" and 
"isotropic" correspond to simulations with globally-averaged temperatures and uniform irradiation, respectively (see text). The dashed curve ("Simulated 
(full)") is taken from a globally-averaged, 1500-day simulation where the first 500 days of initialization have been disregarded, while the dotted curves 
("Simulated") are taken from only 1 day of integration (with no output being disregarded). 



Evaluating equation | |3U in conjunction with equation l |39l l, we get 



where we have defined r = kot. Note that 



|2 + Vs-jexp (-\/37or) + ^ [l - exp (^^Vs-yor^] + ^^[^^ " + Vs-jot^ exp (^-\/37or)] | 



(41) 



70 



is defined such that it is constan t and has no d ependence on either P or e. When e = 
i4\\ reduces to equation (29) of lGuilloll hoid) . 

It follows that the temperature at depth (Tint = K, r 2> 1, 70 ~ 1) becomes 



Tex 



1 

70 



^/3 + 



2(e- 



1, we recover kl = kq, tl = t, 7 



1/4 



KSXo 



(42) 

70 and equation 
(43) 



The additional term in equation ( I43t represents the increase in temperature at depth due to the effect of collision-induced absorption. It is 
non-negligible only when e 3> 1. Specifically, when 



2e 



= rs„ = 1401 



Po 



0.006 cm2 g-i 220 bar/ V9.42 m s 



9p 



(44) 



Even with a wide range of values in e, we get T^o ~ 1500-1800 K. For example, setting e = 10, 100, 1000 and 2000 corresponds to 
Too ~ 1541, 1557, 1699 and 1824 K, respectively. For these respective values of e, collision-induced absorption dominates the long- wave 
optical depth at pressures greater than 24, 2.2, 0.22 and 0.11 bar. Figure |9] plots the temperature-pressure profiles corresponding to these 
values of e, as well as three examples of profiles with different values of 70. We have chosen e — 2000 for the three cases with different 70 
values so as to emphasize the differences between them (which are enhanced with increasing e). 

Note that in deriving the results in th is sub-section, we have retained the values of the first and second (longwave) Eddington coefficients 
(1/3 and 1/2, respectively) as adopted bv lGuillol as well as the closure relation that the ratio of the second to the zeroth moments of 

the sh ortwave intensity is 1/3. Detailed modelling of the opacity sources present in hot lovian atmospheres (e.g.. lGustafsson & Frommhold 
12003") will inform us about the actual value of e, but for now we are content with being able to account for th e effect of collision -induced 
absorption via a single parameter. Readers interested in broader generalizations of iGuilloli(i2010i) are referred to Heng et al. ( 201 Ibh . 



4.4 Baseline models 

4.4.1 Constant opacities (ns = n-L ~ 1) 

As a first step, we test our computational setup against the ns = riL = 1 analytical models of Guillot ( 2010|) . as stated in equations (29) 
and (49) of that study. (See also HansenlEoOSl) . To be consistent with the constant shortwave and longwave opacities of ks — 0.006 cm^ g^^ 
and kl = 0.01 cm g~ , respectively, as adopted bv lGuillol J2OI0I) , we use = 1401 and tlq = 2335 such that 7 = 70 = 0.6. One can 
choose either e = 1 or /; = 0. In the absence of robust empirical constraints, we assume that the longwave optical depth does not depend on 
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Figure 11. Temperature (top panel), zonal wind velocity (middle panel) and meridional wind velocity (bottom panel), averaged over the entire til = 1, hot 
Jupiter simulation domain, as functions of the simulation time (in Earth days). For the velocities, we have computed the root mean square (RMS) values as 
functions of time. The dotted vertical line marks the end of the initialization ("spin up") phase of the simulation. 



latitude. 



(45) 



an assumption we will retain for the rest of the study. 

Figure [Tol compares the analytical versus simulated temperature-pressure profiles. The isotropic profile is obtained from a simulat ion 
where the stellar irradiation is uniform everywhere on the surface. This is equivalent to assuming / = 1 in equation (29) o f iGuillotllioiol) — 
which corresponds to the situation where the exoplanet is uniformly irradiated with the substellar intensity everywhere — and accounts for 
the higher temperature (Too ~ 2129 K) at depth. For the profiles labelled "global", the stellar irradiation is described by equations (|5) and 
l l23t and the resulting temperature-pressure profiles are averaged over the entire globe, thus resulting in a lower temperature (Too ~ 1539 
K) at depth. The temperatures become isothermal at depths at which the stellar iiTadiation becomes completely absorbed (rs > 1). We have 
extracted temperature-pressure profiles only from the first Earth day of integration, such that the zonal winds in the simulations have not had 
time to increase to their quasi-equilibrium values (labelled "Simulated"). Recall that our simulations are initiated with a constant temperature 
Tinit and not from a T-P profile in radiative equilibrium (see i]4.3t . Thus, the reasonable agreement between the dotted and solid curves 
in Figure [Tol demonstrates that the code is implementing the initial radiative condit ion correctly — it is not a rigourous test of the radiative 
transfer scheme. The agreement is imperfect because the analytical expressions of iGuillol 1 201(i) use approximate closure conditions (i.e., 
the Eddington approximations) and also because the vertical grid points in the simula ted profiles are taken to be the larger of the pressure 
half levels in a Simmons-Burridge scheme (see §3.1 o f lHeng. Menou & PhilliDPsl201 l|). The dashed curve in Figure[TO](labelled "Simulated 
(full)") is the globally-averaged temperature-pressure profile where the simulation is executed for 1500 days with the first 500 days being 
disregarded. Its general agreement with the analytical profile constitutes a weak test of the radiative transfer scheme at lower pressures: at 
P < 10 bar, the radiative time scale is t^d < 20 days, which is more than an order of magnitude shorter than our nominal initialization 
period of 500 days, implying that radiative equilibrium has been attained at these pressures. At higher altitudes (P < 10 bar), the fact that 
the simulated profile (dashed curve) in Figure[TO]is generally ~ 50-100 K cooler than the analytical one (solid curve) suggests that thermal 
energy has been converted to mechanical energy in the form of winds, predominantly in the zonal and meridional directions. 

Figure [TT] shows the temperature, zonal wind velocity and meridional wind velocity, averaged over the entire simulation domain, as 
functions of the simulation time. The zonal wind is absent during the first few Earth days of simulation, but gradually builds up as the simu- 
lation attains quasi-equilibrium (sometimes termed "spin up" in the literature). When the zonal wind is accelerated from rest, its equilibrium 
value is reached when the accele rations due to vertical and h orizontal eddy momentum convergence attain steady values, which is expected 
to happen within 100 Earth days I showman & Polvanilbol ih FI It is apparent that 500 days is a reasonable period for the initialization phase. 
The mean temperature plateaus at about 1430 K, which is approximately the equilibrium temperature of the exoplanet. It is possible that the 
very deepest atmospheric layers (P > 100 bar) have not completely "spun up" yet, but in simulations which we executed for longer than 
1500 days we witness that the qualitative trends in our Held-Suarez mean flow quantities are invariant, although the associated quantitative 
details may change slightly. 

An uncertainty of our technique is the value of Cint to adopt: for fccon ~ 0. 1-10'^ W ^ m^ ^ , we haveCint ~ 10^-10^ J K"^ m"^ We 

and 10^ J m~^. The resulting mean temperature-pressure profiles are virtually identical 



thus executed simulations with Cint = 10 , 10 



* The spin up of the zonal wind is shown in Figure 1 1 of lShowman & Polvanil jlOllh , but this specific three-dimensional model was constructed for the case 
of HD 1 89733b. However, we do not expect dramatic changes in the spin up period for the case of HD 209458b. 
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with differences of at most ~ 1 K. We further verified that simulated profiles associated with Cint ~ 10* J 
m~^ deviate substantially (> 50 K) from the analytical ones at depth. 

The zonal-mean zonal wind, temperature, potential te mperature and Euleria n mean streamfunction are shown in Figure [12] It is re- 
assuring that the zonal wind profile resembles Figure 5 of IShowman et al. where a super-rotating, equatorial jet reaches down to 
about 10 bar, flanked by counter-rotating jets at mid-latitudes. The super-rotating and counter-rotating jets have maximum speeds of about 
5.9 km s~^ and -0.6 km s~^, respectively. From examining the potential temperature profiles, one can infer that the atmosphere is baroclinic 
for P < 10 bar (and barotropic at greater pressures), in agreement with the purely dynamical results presented in Figure |8] A difference 
from Figure [8] is that there are now two pairs of circulation cells in each hemisphere (as is evident from examining the Eulerian mean 
streamfunction). We note that convective adjustment is negligibly invoked, which is unsurprising given the deep radiative zone — which 
encompasses both the baroclinic and barotropic atmospheric components — of hot Jupiters. 



4.4.2 Collision-induced longwave absorption (ni^ — 2, e > 1) 

Collision-induced absorptionlfl a phenomenon first discovered by Herzberg ( 1952 ) to be at work within the atmospheres of Neptune 
and Uranus, becomes a non-negligible effect at large pressures. We account for this effect using the formalism derived in ^4.3\ As before, 
our fiducial model assumes tsq ~ 1401, but we now have ro = 2335 and = ero where e > 1. We first check that our simulated and 
hemispherically-averaged (day side only) temperature-pressure profiles are consistent with the analytical ones as derived in ^43] (Figure |9). 
The reasonable agreement between them for e — 10-2000 again demonstrates that our code is implementing the initial radiative condition 
correctly (but does not constitute a rigourous test of the radiative transfer scheme). 

We next focus on a specific value of the correction factor to the longwave optical depth, at the bottom of the simulation domain, due 
to collision-induced absorption: e — 2000. Figure [13] shows the temperature-pressure profiles from our analytical formalism (solid curve), 
a 1-day simulation meant to mimic the case of no atmospheric dynamics (dotted curve) and a 1500-day simulation with the first 500 days 
being disregarded (dashed curve). The simulated profiles are hemispherically-averaged over the day side only. The agreement between the 
solid and dotted curves is reasonable considering that equation ( I41t was formulated for isotropic stellar irradiation (with a dilution factor /). 
We have executed simulations with Cint = 10^, 10^ and 10^ J K^^ m^^ and witnessed differences of at most 1 K in the globally-averaged 
temperature-pressure profiles, again demonstrating that uncertainties in the value of Cint do not significantly affect our results. 

The Held-Suarez mean flow quantities from the full simulation are shown in Figure [T4l where it is apparent that they are qualita- 
ti vely similar to those from Figurell2l The zonal-mean zonal wind profile resembles that obtained from the purely dynamical simulation 
of iHeng. Menou & Phillipp^ Hoi l|) (see the top left panel of their Figure 12). An equatorial, super-rotating jet with a maximum speed of 
about 6.1 km s~^ is again flanked by slower, counter-rotating jets (—0.8 km s^^) at mid-latitudes. At P > 10 bar, the zonal wind structure 
becomes uniform across latitude, which is consistent with the uniform temperature and barotropic structure present. Two pairs of circulation 
cells extend from the equator to the poles on both the day and night sides, in disagreement with the purely dynamical results, presented in 
Figure[8] where only a single pair of cells exists. Unlike in the terrestrial atmosphere, our baseline model shows a baroclinic upper atmosphere 
sitting atop a barotropic lower atmosphere. It is thus somewhat awkward to prescribe the terms "stratosphere" and "troposphere" to the upper 
and lower atmospheres, respectively. The analogue of the tropopause sits at P ~ 10 bar. 

Figure [Is] shows temperature maps as functions of latitude and lo ngitude at P = 2 . 13 mb ar, 216 mbar, 4.69 bar and 21.9 bar. These 
pressure levels are chosen to match those shown in Figure 8 of Heng, Menou & Phillippi ' 2011 ). At P = 2.13 mbar, the location where the 
temperature is at its maximum — w hich we term the "hot spot" — is shifted away from the substellar point, which is so mewhat different 



from the top left panel of Figure 8 of iHeng. Menou & PhillippsI ( 1201 If) and the second panel (from the top) of Figure 16 of ' Showman et al 



(1200 9). Otherwise, Jhe char acteristic chevron-shaped feature, which is seen in all of the three-dimensional simulations of HD 209458b (e.g., 
[siiowman et al.ll2()09i : iRauscher & Menotj|201ol : iHeng. Menou & PhilliDPsir201lh . appears at P = 216 mbar. Deeper in the atmosphere, the 
flow becomes dominated by advection and longitudinal differences in temperature become less apparent. 



4.5 Comparison of different models 

We next compare hemispherically-averaged temperature-pressure profiles from different models in Figure[T6] separated into the day and 
night sides of the simulated hot Jovian atmospheres. We see that for a fixed value of 70 = 0.6, the models with e — 1-1000 essentially have 
the same temperature-pressure profiles, both on the day and night sides, at and above the longwave photosphere. The maximum zonal and 
meridional wind speeds are virtually unchanged: 5.9 km s~^ and from —0.6 to —0.7 km s~^, respectively (at least for the fiducial magnitude 
of the hyperviscosity adopted in this study). The day-night temperature contrasts are almost indiscernible, even though the temperatures 
at depth are modifie d differently due to the ef f ect of collision-induced absorption. The day-night temperature contrasts from the purely 
dynamical model of iHeng. Menou & PhillippsI bOllI) are somewhat different, but we should keep in mind that this simulation employs 



^ For example, the hydrogen molecule in isolation possesses no dipole moment, which implies it cannot absorb photons. However, it may fonn transient 
"super molecules" with other hydrogen molecules under high pressure and the resulting, weak dipole moment allows for absorption. 
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Figure 12. Held-Suarez mean flow quantities for tlie baseline liot Jupiter simulation with nL = 1 (see i|4.4.H . Top left panel: zonal wind (m s ^). Top right 
panel: temperature (K). Middle row: potential temperature (K). Bottom row: Eulerian mean streamfunction (x 10^"^ kg s~^). 



Newtonian relaxation via an equilibrium temperature-pressure profile, which in turn contains a parameter ATeq that sets the (initial) day- 
night temperature contrast as a function of pressure (i.e., ATeq = 530-1000 K from P = 10 bar to P ^ 1 mbar). In our improved 
simulations, the day-night temperature contrasts are computed self-consistently. 

The mod el with 70 = 2 (rsn = 4670, tq = 2 335) exh i bits te mpe rature inversion s on both the day and night sides, consistent with 
the finding bvtaubenv. Burrows & Su darskvl fcoosh . iHansenl (l 2008h and lGuillotI j2O10l) that 70 > 1 models should produce temperature 
inversions. (See Burrows & Orton 2010 for a review of the observational evidence for temperature inversions in hot Jovian atmospheres.) 
Curiously, the temperature at depth is Too ~ 1400 K, which is lower than Teq ~ 1432 K. Thus, this model makes a prediction that infrared 
observations which are able to probe the deep, barotropic layers of a hot Jupiter should infer a equivalent-blackbody temperature which is 
lower than the equilibrium temperature of the exoplanet. For the model with 70 = 2, we note that the Held-Suarez mean flow quantities are 
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Figure 13. Compaiing analytical and simulated temperature-pressure profiles for the case of riL = 2 and e = 2000. The solid curve is based on our 
generalization of the Guillot (2010) analytical profile (see fi43\ . The dotted curve ("Simulated") is taken from a hemispherically-averaged. 1-day simulation, 
while the dashed curve ("Simulated (full)") is taken from a hemispherically-averaged, 1500-day simulation where the first 500 days of initialization have been 
disregarded. Note that both of the simulated profiles are extracted from the day side only. 



qualitatively similar to the baseline hot Jupiter models previously presented (except for the temperature field), despite the differences present 
in the temperature-pressure profiles. 

These differences in the temp erature-pressure profiles a re relevant to observations, because the longwave photosphere is, in the absence 
of clouds, located at a pressure of JSurrows & Liebertlll993l) 

^ ^ ^ 63 mbar ( ( " . (46) 

We conclude that the quantity 70 controls the presence or absence of a temperature inversion, as well as the day-night temperature 
contrast. The key differences from the one- or two-dimensional models is that the redistribution of heat/energy and the depths at which 
shortwave/longwave absoiption occurs are calculated self-consistently once the optical depths are specified. We also note that while we 
have shown hemispherically-averaged temperature-pressure profiles for the day and night sides, our calculations generally produce three- 
dimensional T-P profiles. 



5 DISCUSSION 
5.1 Observational consequences 

A direct consequence of our results concerns the vertical mixing of atmospheric fluid. As explained in equat ion ij2.6l the strength o f 



atmospheric circulation is quantified by the Eulerian mean streamfunction. On Jupiter, we have * ~ 10"' kg s"^ \Uu & Schneideilboiol) . 
As shown by our results in Fi gures |8] 1 1 2 1 and [T4l the circulation in hot Jovian atmospheres has a strength of ^l/ ~ 10^* kg s~^, which is 
~ 10'* times stronger than on Jupiter. Since the circulation cells extend from ~ 1 mbar to ~ 10 bar in our models, the relevant length scale is 
~ lOJy, where H = k-^T /fhgp is the pressure scale height. If we pretend that we may prescribe a "diffusion coefficient" to the atmospheric 
circulation, then 

7^..~10//y^^«.($,ao)dao~-^~10 cm s isqq k io14 i ,g ) [2^ 9.44 x 10^ k m 220 bar J ' ^^'^ 
which is consistent with the required value of Kzz ~ 10^-10^^ cm^ s^^, estimated by Spiegel. Silverio & Burrows I dlOOgh for HD 209458b, 



in order to keep particles associated with TiO, with sizes of 0.1-10 /im, aloft such that they may act as shortwave absorbers capable of pro- 
ducing a temperature inversion. Unfortunately, the Eulerian mean streamfunction is not a good representation of the vertica l or horizontal 
mixing of particles embedded in the flow (termed "tracers" in the atmospheric science community; see §12 of VallisI 200^ . So while the 



enhanced vertical mixing (relative to Jupiter), via large-scale circulation cells present in hot Jovian atmospheres, may be a tempting and 
plausible mechanism f or TiO to be maintained a t high altitudes, a more detailed study of the interaction of tracers with the atmospheric circu- 
lation is required. (See Youdin & Mitchelll201(]| for toy models of vertical, turbulent mixing in hot Jovian atmospheres and their implications 



for the exoplanetary radii.) 

An obs ervational way of distinguishing between different hot Jovian atmospheres is to examine their thermal phase curves jCowan & Agol 



200 sL 2011 ). The flux associated with the outgoing longwave radiation (OLR), denoted by J-qlr = Jx)lr(Q, is the emergent flux from 
the longwave photosphere. In Figure [17] we show maps of J-qlr for two models (til — 2, e — 2000; 70 = 0.6 versus 2), where it is 
apparent that the chevron- shaped feature seen at P ~ 0.1 bar resides near the longwave photosphere. The OLR flux may be used to construct 
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Figure 14. Same as Figure[T2] but for a simulation with nL = 2 and e = 2000. 



the thermal phase curves, as was done for exo-Earths by ISelsis. Wordsworth & Forgeli ( 120111) . One has to first average the OLR flux over 
latitude. 



1 

TrJ-^ 

The flux associated with the thermal phase curve, J-'phasc, can then be constructed using equation (7) of ICowan & Ago 



(48) 



(^olr) = Aq + Aj cos + Bj sin (j» , 



(49) 



.T^phasc = 2Ao + ^ Cj COS (ji)) + Dj sin (j't/j) , 
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Figure 15. Snapshots of the temperature field as functions of latitude ($) and longitude (©). for the hot Jupiter simulation with riL = 2 and e = 2000, at 
P = 2.13 mbar (top left panel), 216 mbai' (top right panel), 4.69 bar (bottom left panel) and 21.9 bar (bottom right panel). The snapshots are taken at 1500 
days after the start of the simulation. Temperatures are given in K. 
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Figure 16. Hemispherically-averaged temperature-pressure profiles, separated into the day and night sides, for the hot Jupiter models presented in this study. 
Left panel: comparing rtL = 1 and e = 10, 100 and 1000 cases. Note that the rtL = 1 (e = 1) and e = 10 cases are almost indistinguishable in this plot. 
Right panel: comparing profiles with e = 2000 and different values of 70. The curves labelled by "HMP" are taken from the HD 209458b model of Heng, 
Menou & Phillipps (201 1), which only considers atmospheric dynamics and uses Newtonian relaxation to mimic radiative cooling. The dashed, horizontal line 
indicates the approximate location of the longwave photosphere. In the right panel, the condensation curves for TiO from Spiegel, Silverio & Burrows (2009) 
are shown. 
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Figure 18. Latitudinally-averaged outgoing longwave radiation (OLR; top panel) as a function of the longitude © and thermal phase curve (bottom panel) as a 
function of the phase angle tj), for the two models shown in Figure [Tt] In the bottom panel, the thermal phase curves peak before the secondary eclipses occur 



where = + vr, the phase angle is denoted by i/) and the summation coefficients are related via the following expressions: 



2 



and 




(50) 

i = i, 

(51) 

J = l- 

We perform the fits, using the function described in the first expression in equation ( I49t . to our computed { J-olr) curves by considering terms 
up to n = 13. We then obtain J-"Dhase by transforming the Aj and Bj coefficients, using equations l ISQt and | |5U . to Cj and Dj coefficients. 
As pointed out bv lCowan & Agol 1 20081) . the odd sinusoidal modes do not contribute to the thermal phase curve: we set Cj = Dj — when 
j is odd except for j = 1. This implies that information is lost when the latitudinally-averaged OLR flux is converted to the thermal phase 
curve. For example, for the 70 = 0.6 model, we have [{A2/Aof + iB2/Aof]^^^ ~ 0.08 while [{As/Aof + {Ba/AofY^^ « 0.04. For 
the 70 = 2 model, we have [{A2/Aof + {B2/AofY^^ « 0.15 while [{Aa/Aof + [Ba/Aof]^^^ « 0.09. These estimates suggest that 
the second lowest odd mode (j — 3) has a small but non-negligible contribution to the OLR map. It is worth noting that the odd modes 
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Figure 19. Same as Figure [78] but for comparison of the 70 = 0.6 and t = 2000 model with hyperviscous time scales of f ^ = 10 ^,10 ^ and 10 ^ HD 
209458b day. The hot spot offset appears to be somewhat robust to our ignorance of ti,. 



do not contribute to the thermal phase curve only when the following assumptions are valid: the exoplanet rotates edge-on; there is no limb 
darkening associated with the exoplanet; and there is no temporal variability associated with the OLR. 

In FigurefTSl we calculate (J-olr) and J-'phaso for the 70 = 0.6 and 70 = 2 models previously shown in Figure[T7] For the latitudinally- 
averaged OLR flux (top panel), the angular offs ets of the peaks from th e substellar point (O = 180°) are 33.8° and 26.3° for the 70 — 0.6 
and 70 = 2 models, respectively. As realized bv lCowan & Agoll (1201 in . the corresponding offsets in the thermal phase curve (bottom panel) 
are somewhat larger: 39.4° and 31.9°, which correspond to temporal offsets of 9.3 and 7.5 hours, respectively, easily discernible with current 
technology. Notice that the phase angle 1/) is defined in such a manner that the peaks of the thermal phase curves occur before superior 
conjunction (i.e., secondary eclipse). 



We may be concerned that drag within the atmospheres may alter the location of the hot spot, as implied by the study of lShowman & Polvani 



1 20111) . Within the con text of the simulations, drag may manifest itself in two forms: physical (e.g., Ohmic) versus numerical (i.e., horizontal 
dissipation; see §3.3 of^ Heng. Menou & PhillipDsl201 1 ). We restrict our discussion to the uncertainties associated with numerical drag, which 
takes the form of hyperviscosity in our spectral simulations. To illustrate the effects of numerical drag, we execute two more simulations 
with 70 — 0.6 and e — 2000 (our default model), but with the time scale associated with hyperviscosity (t^) decreased by factors of 10 
and 100. The corresponding latit udinally-averaged OLR flux and thermal phase curves are shown in Figure [19] As already demonstrated in 



irresponoin g 
: Phillipp j (2 



Heng. Menou & Phillippa (20 HI) , changing the hyperviscosity alters the zonal wind speeds: the maximum and minimum speeds of 6.1 km 



s ^ and —0.8 km s ^, respectively, from our default model now become 5.8 km s"^ and -1.2 km s"^ for U = 10'^ HD 209458b day 
and 5.5 km s~^ and —1.1 km s~^ for ti, — 10~^ HD 209458b day. Since hyperviscosity is a numerical tool and it is difficult to specify 
its magnitude from first principles, zonal wind speeds in hot Jovian atmospheres are not robust predictions of general circulation models, 
at least at the < 10% level. However, the hot spot offset appears to be a somewhat robust prediction of the models: for — , 
and 10~^ HD 209458b day, the temporal offsets are 9.3, 9.7 and 9.7 hours, respectively. Therefore, despite our ignorance of the hot spot 
offsets measured from thermal phase curves offer an opportunity to distinguish between hot Jovian atmospheres characterized by different 
values of 70 . 



5.2 Summary 

The salient points of our study can be summarized as follows: 

• We have successfully implemente d a simple (dual-band, two-stream ) radiative transfer scheme in conjunction with the FMS spectral 
dynamical core previously described in iHeng. Menou & Phillipp^ 1 2011 ). We have tested our computational setup by simulating the atmo- 
spheres of Earth-like (exo)planets and hot Jupiters. Our models possess an intermediate degree of sophistication in a required hierarchy of 
three-dimensional models of atmospheric circulation. 

• We have generalized the analytical formalism of lGuillot to calculate temperature-pressure profiles of hot Jovian atmospheres, in 
radiative equilibrium, which include the effect of collision-induced absorption via a single parameter e. 
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• Atmospheric circulation on hot Jupiters is characterized by the presence of large circulation cells extending from the equator to the 
poles, due to the effect of Coriolis deflection being weakened by slower rotation. While a transitional layer ("tropopause") may exist between 
the baroclinic and barotropic components of the atmosphere, the former is generally located at greater altitudes while the latter resides deeper 
in the atmosphere. This is the opposite from the situation on Earth (and a hypothetical, tidally-locked Earth), where the stratosphere sits on 
top of the troposphere. 

• Large-scale circulation cells in hot Jovian atmospheres offer a plausible mechanism for maintaining TiO at high altitudes, such that it 
may act as a shortwave absorber capable of producing a temperature inversion. However, a more detailed study of the interaction between 
particles embedded in the flow ("tracers") and the atmospheric circulation is required before any robust conclusions can be drawn. 

• The absence or presence of a temperature inversion in the baroclinic component of a hot Jovian atmosphere, as well as the temperature 
contrast between the day versus night sides, is determined by the ratio of shortwave to longwave opacity normalizations (70), at least in a 
cloud- or haze-free scenario. 

• The angular offset between the hot spot — associated with the equatorial, super-rotating wind expected to be present in hot Jovian 
atmospheres — and the substellar point may help us distinguish between hot Jovian atmospheres characterized by different values of 70, and 
appears to be robust to our ignorance of hyperviscosity. 



6 ACKNOWLEDGMENTS 

K.H. acknowledges support from the Zwicky Prize Fellowship and the Star and Planet Formation Group at ETH Zurich. D.M.W.F. 
acknowledges support from NSF grants ATM-0846641 and ATM-0936059. K.H. thanks Kristen Menou, Frederic Pont, Nikku Madhusud- 
han, Hans Martin Schmid, Tim Merlis, Eric Agol, David Sing, Dave Spiegel, Michael Meyer, Emily Rauscher, Tapio Schneider and Paul 
O'Gorman for useful conversations. We additionally thank Nick Cowan, Andrew Youdin and Rene Heller for useful comments following 
their reading of an earlier version of the manuscript. K.H. is especially grateful to Kristen Menou for providing invaluable guidance during 
the initial stages of this study and Nick Cowan for providing a crash course on thermal phase curves. We thank the anonymous referee for 
constructive criticism which improved the quality of the manuscript. This work benefited from competent and dependable technical support 
by Olivier Byrde, Eric Miiller, Adrian Ulrich et al. of the Brutus computing cluster team, as well as steadfast administrative support by 
Marianne Chiesi of the Institute for Astronomy, at ETH Zurich. 



REFERENCES 

Andreas, E.L., Claffey, K.J., Jordan, R.E., Fairall, C.W., Guest, PS., Persson, P.O.G., & Grachev, A.A. 2006, Journal of Fluid Mechanics, 
559, 117 

Berger, A., & Loutre, M.F. 1991, Quaternary Science Reviews, 10, 297 

Betts, A.K., 1986, Quarterly Journal of the Royal Meteorological Society, 112, 677 

Betts, A.K., & Miller, M.J. 1986, Quarteriy Journal of the Royal Meteorological Society, 112, 693 

Brown, T.M., Charbonneau, D., Gilliland, R.L., Noyes, R.W., &. Buitows, A. 2001, ApJ, 552, 699 

Burrows, A., & Liebert, J. 1993, Reviews of Modem Physics, 65, 301 

Burrows, A., & Orton, G. 2010, Giant Planet Atmospheres, in "Exoplanets", ed. S. Seager, p. 419^40 (Tucson: University of Arizona 
Press) 

Cooper, C.S., & Showman, A.P 2005, ApJ, 629, L45 
Cooper, C.S., & Showman, A.P 2006, ApJ, 649, 1048 
Cowan, N.B., & Agol, E. 2008, ApJ, 678, L129 
Cowan, N.B., & Agol, E. 2011, ApJ, 726, 82 
Crowley, T.J. 2000, Science, 289, 270 

Dobbs-Dixon, I., Gumming, A., & Lin, D.N.C. 2010, ApJ, 710, 1395 

Frierson, D.M.W., Held, I.M., & Zurita-Gotor, P. 2006, Journal of the Atmospheric Sciences, 63, 2548 

Frierson, D.M.W., Held, I.M., & Zurita-Gotor, P. 2007, Journal of the Atmospheric Sciences, 64, 1680 

Frierson, D.M.W., Lu, J., & Chen, G. 2007, Geophysical Research Letters, 34, L18804 

Frierson, D.M.W. 2007a, Journal of the Atmospheric Sciences, 64, 1959 

Frierson, D.M.W. 2007b, Journal of the Atmospheric Sciences, 64, 2076 

Frierson, D.M.W. 2008, Journal of the Atmospheric Sciences, 65, 1049 

Gartatt, J.R. 1994, Earth-Science Reviews, 37, 89 

Gordon, C.T., & Stem, W.F. 1982, Monthly Weather Review, 110, 625 

Guillot, T. 2010, A&A, 520, A27 

Gustafsson, M., & Frommhold, L. 2003, A&A, 400, 1161 
Hansen, B.M.S. 2008, ApJS, 179, 484 

Held, I.M., & Suarez, M.J. 1994, Bulletin of the American Meteorological Society, 75, 1825 



© 2013 RAS, MNRAS 0Q0.[Tll28] 



Atmospheric circulation of exoplanets II 27 



Held, I.M. 2005, Bulletin of the American Meteorological Society, 86, 1609 
Heng, K., Menou, K., & Phillipps, P.J. 2011, MNRAS, 413, 2380 
Heng, K., & Vogt, S.S. 2011, MNRAS, 415, 2145 

Heng, K., Hayek, W., Pont, R, & Sing, D.K. 201 1, preprint i arXiv:1107.1390l > 
Herzberg, G. 1952, ApJ, 115, 337 

Holton, J.R. 2004, An Introduction to Dynamic Meteorology, 4th edition (Massachusetts: Elsevier) 
Hubbard, W.B. 1968, ApJ, 152, 745 

Hubbard, W.B., Burrows, A., & Lunine, J.I. 2002, Annual Reviews of Astronomy & Astrophysics, 40, 103 
Hubeny, I., Burrows, A., & Sudarsky, D. 2003, ApJ, 594, 1011 
Iro, N., Bezard, B., & Guillot, T. 2005, A&A, 436, 719 

Kang, S.M., Held, I.M., Frierson, D.M.W., & Zhao, M. 2008, Journal of Climate, 21, 3521 
Koskinen, XT, Aylward, D., Smith, C.G.A., & Miller, S. 2007, ApJ, 661, 515 
Langton, J., & Laughlin, G. 2008, ApJ, 674, 1106 

Lewis, N.K., Showman, A.P, Fortney, J.J., Mariey, M.S., Freedman, R.S., & Lodders, K. 2010, ApJ, 720, 344 
Liu, J., & Schneider, T. 2010, Journal of the Atmospheric Sciences, 67, 3652 
Lorenz, E.N. 1955, Tellus, 7, 157 

Lu, J., Chen, G., & Frierson, D.M.W. 2010, Journal of the Atmospheric Sciences, 67, 3984 
Manabe, S., Smagorinsky, J., & Strickler, R.F. 1965, Monthly Weather Review, 93, 769 
Mazeh, T., et al. 2000, ApJ, 532, L55 

Merits, T.M., & Schneider, T. 2010, Journal of Advances in Modeling Earth Systems — Discussion (JAMES-D), in press 

Mihalas, D. 1978, Stellar Atmospheres, 2nd edition (San Francisco: Freeman) 

Mitchell, J.L., Pierrehumbert, R.T., Frierson, D.M.W., & Caballero, R. 2009, Icarus, 203, 250 

O'Goi-man, PA., & Schneider, T. 2008, Journal of Climate, 21, 3815 

Palluconi, RD., & Kieffer, H.H. 1981, Icarus, 45, 415 

Pauluis, O., Czaja, A., & Korty, R. 2008, Science, 321, 1075 

Peixoto, J.P., & Oort, A.H. 1984, Reviews of Modem Physics, 56, 365 

Phillips, N.A. 1957, Journal of Atmospheric Sciences, 14, 184 

Pierrehumbert, R.T., & Swanson, K.L. 1995, Annual Reviews of Fluid Mechanics, 27, 419 
Pierrehumbert, R.T. 2010, Principles of Planetary Climate (New York: Cambridge University Press) 
Ramanathan, V., & Coakley, Jr., J.A. 1978, Reviews of Geophysics and Space Physics, 16, 465 
Rauscher, E., & Menou, K. 2010, ApJ, 714, 1334 

Selsis, F, Wordsworth, R.D., & Forget, F. 2011, A&A, in press ( arXiv:lI04.4763V 2) 
Showman, A.P, & Guillot, T. 2002, A&A, 385, 166 

Showman, A.P, Cooper, C.S., Fortney, J.J., & Mariey, M.S. 2008, ApJ, 682, 559 

Showman, A.P, Fortney, J.J., Lian, Y., Mariey, M.S., Freedman, R.S., Knutson, H.A., & Charbonneau, D. 2009, ApJ, 699, 564 

Showman, A.P, & Polvani, L.M. 2010, Geophysical Research Letters, 37, L1881 1 

Showman, A.P, & Polvani, L.M. 201 1, preprint (arXiv: 1 103.3101i'I) 

Spiegel, D.S., Silverio, K., & Burrows, A. 2009, ApJ, 699, 1487 

Stevenson, D.J. 1982, Annual Reviews of Earth and Planetary Science, 10, 257 

Thi-astai-son, H.Th., & Cho, J.Y-K. 2010, ApJ, 716, 144 

Troen, I.B., & Mahrt, L. 1986, Boundary-Layer Meteorology, 37, 129 

Vallis, G.K. 2006, Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation (New York: Cambridge Univer- 
sity Press) 

Washington, W.M., & Parkinson, C.L. 2005, An Introduction to Three-Dimensional Climate Modeling, second edition (Sausalito: Univer- 
sity Science Books) 
Xu, K.-M., & Emanuel, K.A. 1989, Monthly Weather Review, 117, 1471 
Youdin, A., & Mitchell, J.L. 2010, ApJ, 721, 1 113 



APPENDIX A: ATMOSPHERIC BOUNDARY LAYER 

There are two main pieces of key physics which need to be included within the atmospheric boundary layer. Firstly, the turbulence 
generated within the layer acts as a form of drag between the surface and the atmosphere. Secondly, since the turbulent eddies generated are 
expected to be three-dimensional, energy flows towards the smallest length scales (i.e., a forward Kolmogorov cascade) where it is dissipated 
by molecular viscosity. The process of turbulent mixing may be described by the diffusion equation. Therefore, a boundary layer scheme 
needs to approximately capture the effects of drag and diffusion. 

Operationally, drag between the surface and the lowest model layer is expressed in the form of drag coefficients acting on the temperature 
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and velocit y fields, which are calculated using th e "Monin-Obukhov" drag laws, considered standard in the atmospheric science community 
(see §2c of Frierson. Held & Zurita-Gotoij|2006 ): 



1 - 



Tlzo < 0, 

< TZif, < 7?.i,crit, 
crit 5 



where 



C = ln 



zo 



-rough 



(Al) 



(A2) 



and Zo is the height of the lowest model layer. It should be noted that equation iAl\ is considered a simplified Monin-Obukhov formulation 
— more sophisticated functions exist to describe Cmo- The key physics is captured in three numbers: trough, TZi and KvK. The first of these 
is the roughness length trough, which parametrizes the macroscopic effects of the surface type on the drag. Within the Monin-Obukhov 
framework, it is the vertical height at which t he horizontal win ds go to zero. For example, it is trough ~ 10~* m over open water and ~ 1 
m over urban terrain. Following lFrierson. Held~& Zurita-Gotoj 120061) , we adopt trough = 3.21 x 10~^ m for our baseline Earth-like model 
( CT . 



Denoting the horizontal wind speed by v, the second of these numbers is the bulk Richardson number (e.g., §3.7.1 of Washington & ParkinsonI 

2005h , 
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5p ddT ( dv 
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which quantifies the effects of shear versus vertical stability. When the bulk Richardson number falls below a critical value T^^.crit ~ 1-10, 
the atm osphere is turbulent. Above this critical value, the drag coefficients are set to zero. Again following iFrierson. Held & Zurita-Gotoi 
(200i), we set 7?.i,crit ~ 1 for our Earth-like model. In equation ( lAU . the quantity TZig is the bulk Richardson number evaluated at the lowest 
model layer. 

The third of these number s is the von Karman constant, which from laboratory experiments and terrestrial atmospheric observations is 
estimated to be Kvk ~ 0.4 (see Andreas et al.ll2006l and references therein). We thus adopt Kvk ~ 0.4 in our Earth-like simulations. 



T o treat turbulent diffusion, we first determine the boun dary layer depth fe by c alculating the height where TZi = 7?.i,crit JFrierson. Held & Zurita-Gotor 
200d) F°lThe turbulent diffusion coefficient is computed as fcoen & Mahr3ll986l) 



/C 




z < z 



Btt) z'<z<h, 



(A4) 



where 



z' = fth. (A5) 

The boundary layer fraction defines a "surface layer" of height z' which is assumed to be in heat and momentum equilibrium with the 
surface — this formulation is made in the intere st of computational efficiency, so as to avoid an iterative calculation to establish equilibrium. 
Following lFrierson. Held & Zurita- Gotorl(l2006l) . we adopt ft = 0.1. Above the surface layer, the diffusion coefficient smoothly goes to zero 
at the top of the boundary layer (z = h). The surface layer diffusion coefficient is 



K.{z) 



^1/2 

KvK VQZ C,;,-, 



Uj ln(z/z„„gh) 



7^^o < 0, 

7^,o > 0, 



(A6) 



where vq represents the horizontal velocity in the lowest model layer. The diffusion coefficients act on the velocity field and the dry static 
(specific) energy, the later of which is defined as 

fdry = CpT + QpZ. (A7) 

The Monin-Obukhov framework describes mixing induced via both shear and convection, and allows for a smooth transition between 
the two regimes. A key difference between it and the convective adjustment scheme ( i]2.4l ( is that the former mixes momentum while the latter 
does not. An obstacle to applying the Monin-Obukhov scheme to the atmospheres of hot Jupiters is that the four parameters are considered 
to be free. 



" At the top of the boundary layer, there often exists a stable layer where the turbulent motions from beneath do not penetrate into, termed the "capping 
inversion". 
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